'PREPRINT       Dr T Clarke W3IWI
'Converted to the Amiga by T Weatherley G3WDI
'Added to by A Hewat G8NTH
'Version 1.00    First Amiga version
'Last mod 26/09/88 by ACH

DIM T$(20),S$(40),I$(40),C(3,2):P=0
DEFDBL A-Z

C8$=CHR$(10)+CHR$(10)+CHR$(10)+CHR$(10)
50 C9$=CHR$(12)+CHR$(7):PRINT C9$
PRINT "AMSAT ORBITAL PREDICTION PROGRAM"
PRINT  "COPYRIGHT 1980  DR. TOM CLARKE W3IWI"
PRINT  "TRANSLATED FOR AMSAT UK BY G3WDI"
PRINT "COPYRIGHT AMSAT-UK. 1988"
REM SET STARTING DATA
INPUT "start year = ",Y:Y=Y/100:Y=INT(100* (Y-INT(Y)) + 0.1)
IF Y/4=INT(Y/4) THEN F9=1 :ELSE F9=0
INPUT "     month = ",M
INPUT "       day = ",D
 DEF FNT$(D)= CHR$(48+INT(D/10))+CHR$(48+D-10*INT(D/10))
T$=FNT$(D)+"/"+FNT$(M)+"/"+FNT$(Y)+" at "
RESTORE daysofyear:
FOR I= 1 TO M:READ D9:NEXT:D8=D+D9
IF M>2 THEN D8=D8+F9 :PRINT "    day = ",D8
PRINT
INPUT"Start : UTC hours= ",H
INPUT"             Mins= ",M
T7=D8+H/24+M/1440
T$=T$+FNT$(H)+FNT$(M)+":00"
PRINT ''
INPUT"duration hours= ",H1
INPUT"          mins= ",M1:T8=T7+H1/24+M1/1440
 PRINT ''
INPUT "TIMESTEP: Min= ",M2:T9=M2/1440
PRINT '
PRINT "    from ",:PRINT  USING"######.####";T7,:PRINT ," to ";:PRINT  USING"######.####";T8
RESTORE Physicaldata:
READ P1,C,R0,F,G0,G1
P2=2*P1:P0=P1/180:F=1/F
Getyear:
READ Y1,G2:IF Y1=Y THEN statdata:
IF Y1>0 THEN Getyear:
IF Y1=0 THEN PRINT "UNABLE TO FIND YEAR:STOP
statdata:
RESTORE stationdata:
READ L9,W9,H9,C$,E8
GOSUB observerdata: PRINT''
DEFDBL Q,K
GOSUB pobsdata:
PRINT "-------------satellite selection menu -------------------------"
I=0:RESTORE keps:
mnu:
READ S$:IF S$="END" THEN choice: :ELSE READ I$:I=I+1
PRINT "Entry # ";I;" for ";S$:PRINT "   ID= ";I$
READ D,D,D,D,D,D,D,D,D,D,D,D,D,D : GOTO mnu:
choice:
INPUT "Select entry # ",J:RESTORE keps:
IF J<1 OR J>I THEN choice:
FOR I= 1 TO J:READ S$,I$,Y3,D3,H3,M3,S3
READ K0,M0,N0,A0,I0,E0,W0,O0,F1:NEXT
PRINT :PRINT  "Doppler calculated for frequency = ",F1," MHz "
INPUT "Change frequency to (0 for default) ",D
INPUT "printer Y/y / n   ",P$ :IF P$="Y" OR P$ ="y" THEN printer=1
IF D<>0 THEN F1=D
IF Y3=Y THEN 500 :ELSE PRINT "ELEMENTS NOT FROM CURRENT YEAR":STOP
500 IF D3=INT(D3) THEN T0=D3+H3/24+M3/1440+S3/86400 :ELSE T0=D3
PRINT :PRINT "Orbital Elements for "'S$
PRINT "    Reference Id.= ",I$:PRINT 
PRINT "Reference Epoch = ";Y3;" + ";T0
PRINT "Starting  Epoch = ";Y;" + ";T7;" = ";T$
PRINT 
PRINT "Parameter";SPC(9);"Reference";SPC(22);"Starting"
T=T7
GOSUB matrix: : GOSUB meananomoly: 
PRINT "Orbit number ";SPC(10);K0;SPC(10);K
PRINT "Mean Anomaly ";SPC(10);M0;SPC(10);M/P0
PRINT "Inclination  ";SPC(10);I0
PRINT "Eccentricity ";SPC(10);E0
PRINT "Mean Motion  ";SPC(10);N0
PRINT "S.M.A.,Km    ";SPC(10);A0
PRINT "Arg.Perigee  ";SPC(10);W0;SPC(10);W
PRINT " R. A. A. N. ";SPC(10);O0;SPC(10);O
PRINT "Freq. Mhz.   ";SPC(10);F1

K9=5.47102e-315:K8=5.47102e-315:

maincomputationloop:
FOR T=T7 TO T8 STEP T9:K7= INT(T):GOSUB meananomoly
IF K=K9 THEN 720:GOSUB matrix::K8=5.47102e-315:K9=5.47102e-315
720 GOSUB Keplerseq: 
GOSUB parameters: 
IF E9<E8 THEN 830 
IF K7=K8 AND K9=K THEN 770
730 IF K=K9 THEN 760 :ELSE GOSUB pager:
K9=K
IF printer=1 THEN 700 :ELSE 740
700 CLOSE #2:OPEN "LPT1:" FOR OUTPUT AS #2
701 PRINT #2," U.T.C.    AZ    EL   DOPPLER  RANGE    HEIGHT   LAT   LONG   PHASE"
702 PRINT #2,"HHMM:SS   deg   deg      Hz      km       km     deg    deg    (256)"

740 'PRINT " U.T.C.    AZ    EL   DOPPLER  RANGE    HEIGHT   LAT   LONG   PHASE"
750 'PRINT "HHMM:SS   deg   deg      Hz      km         km   deg    deg   (256)"
760 PRINT SPC(10);"---DAY# ";K7;"---ORBIT # ";K;"----"
770 K8=K7:T4=T-K7:S4=INT(T4*86400 +0.5):H4=INT(S4/3600+1e-06)
780 M4=INT((S4-H4*3600)/60+1e-06):S4=INT(S4-3600*H4-60*M4)
790 T$=FNT$(H4)+FNT$(M4)+":"+FNT$(S4)
800 F9=-F1*1e+06*R8/C  : DEF FNI(D)=SGN(D)*INT(ABS(D)+0.5)
810 PRINT #2, T$;:PRINT #2,USING "######";FNI(A9);FNI(E9);:PRINT#2, USING"#######";FNI(F9);
PRINT #2,USING "#########";FNI(R5);FNI(R-R0);:PRINT #2, USING"######";FNI(L5);FNI(W5);M9
830 NEXT:PRINT : INPUT"PRESS RETURN TO CONTINUE:",D$:GOTO 50
pager:
PRINT  C9$;C$;"  lat=";L9;"  W.Long=";W9;"  Ht. = "H9;
P=P+1:PRINT SPC(10);" Pg.="P
PRINT  SPC(10);"---MINIMUM ELEVATION= ";E8;"Deg.--"
PRINT  :RETURN


STOP






observerdata:
L8=L9*P0
S9=SIN(L8)
C9=COS(L8)
S8=SIN(-W9*P0):C8=COS(W9*P0)
R9=R0*(1-(F/2)+(F/2)*COS(2*L8))+H9/1000
L8=ATN((1-F)^2*S9/C9):Z9=R9*SIN(L8)
X9=R9*COS(L8)*C8:Y9=R9*COS(L8)*S8:RETURN
pobsdata:
PRINT C$;"   Lat.="L9;" W.Long.=";W9;"   ht.=";H9
PRINT TAB(5)"---Minimum elevation -";E8;"Deg--"
PRINT 
RETURN
matrix:
IF N0>0.1 THEN A0=((G0/(N0^2))^(1/3))
IF N0<=0.1 THEN N0=SQR(G0/(A0^3))
E2=1-E0^2:E1=SQR(E2):Q0=M0/360+K0
K2=9.95*((R0/A0)^3.5)/(E2^2)
S1=SIN(I0*P0):C1=COS(I0*P0)
O=O0-(T-T0)*K2*C1
S0=SIN(O*P0)
 C0=COS(O*P0)
W=W0+(T-T0)*K2*(2.5 * (C1^2)-0.5)
S2=SIN(W*P0)  
C2=COS(W*P0)
C(1,1)=+(C2*C0)-(S2*S0*C1) 
C(1,2)=-(S2*C0)-(C2*S0*C1)
C(2,1)=+(C2*S0)+(S2*C0*C1) 
C(2,2)=-(S2*S0)+(C2*C0*C1)
C(3,1)=+(S2*S1)           
C(3,2)=+(C2*S1) : RETURN
meananomoly:
Q=N0*(T-T0)+Q0:K=INT(Q):M9=INT((Q-K)*256):M=(Q-K)*P2:RETURN


keps:
DATA "Cosmos 1766"
DATA "LES"
DATA 87,230.74062106,0,0,0
DATA 999
DATA 88.5386
DATA 13.8338111720,0
DATA 82.5567
DATA 0.0011
DATA 271.1364
DATA 29.0651
DATA 137.4

DATA "MIR"
DATA "NASA"
DATA 87,350.92010311,0,0,0
DATA 10492
DATA 141.9009
DATA 15.741999605,0
DATA 51.63
DATA 0.0019
DATA 218.0809
DATA 291.7837
DATA 143.625

DATA "METEOR 2/16"
DATA "NASA"
DATA 87,349.94089450,0,0,0
DATA 1657
DATA 58.1840
DATA 13.83330619,0
DATA 82.5557
DATA 0.0011
DATA 301.8209
DATA 294.6429
DATA 137.850

DATA "METEOR 2/15"
DATA "NASA"
DATA 87,349.92970702,0,0,0
DATA 7849
DATA 250.9101
DATA 13.83767605,0
DATA 82.5363
DATA 0.0015
DATA 109.3764
DATA 15.1388
DATA 137.85

DATA"NOAA-10"
DATA"NASA"
DATA 87,346.71759361,0,0,0
DATA 6412
DATA 308.0062
DATA 14.22534037,0
DATA 98.7008
DATA 0.0015
DATA 52.2462
DATA 15.1388
DATA 137.5

DATA "NOAA-9"
DATA "NASA"
DATA 87,341.20109492,0,0,0
DATA 15372
DATA 280.3594
DATA 14.1153643,0
DATA 99.0755
DATA 0.0016
DATA 79.9437
DATA 308.5130
DATA 137.62

DATA "METEOR 30"
DATA "NASA"
DATA 87,346.29912200,0,0,0
DATA 40878
DATA 72.8717
DATA 14.98493528,0
DATA 97.7216
DATA 0.0041
DATA 286.7977
DATA 72.4335
DATA 137.010

DATA "SALYUT 7"
DATA "NASA"
DATA 87,350.91797151,0,0,0
DATA 32398
DATA 207.1197
DATA 15.31746328,0
DATA 51.6122
DATA 0.0001
DATA 152.9965
DATA 94.1426
DATA 20.005
DATA "END"                     
'ADD NEW DATA HERE EPOCH ORBIT MA MM INC ECC ARG RAAN FQ




Physicaldata:
DATA 3.1415926535,2.997925E5,6378.160,298.25:'PI/VEL LIGHT/RAD EARTH/FLATTENING
DATA 7.5369793E13,1.0027379093:'GM OF EARTH/SID/SOLAR TIME
siderealdata:
DATA 86 ,.27
DATA 87,0.2753552680
DATA 88,0.2746921630
DATA 89,0.2767669670
DATA 90,0.2767669760
DATA 0,0
stationdata:
'put your data here lat ,long, height ,call,minimum elevation
DATA 51.26833,0.5633,35,"G8NTH",1
 
 


daysofyear:
DATA 0,31,59,90,120,151,181,212,243,273,304,334,365

Keplerseq:
E=M+E0*SIN(M)+0.5*(E0^2)*SIN(2*M)
10420 S3=SIN(E):C3=COS(E):R3=1-E0*C3:M1=E-E0*S3
M5=M1-M:IF ABS(M5)<1e-06 THEN 10450 :ELSE E=E-M5/R3:GOTO 10420
10450 X0=A0*(C3-E0):Y0=A0*E1*S3:R=A0*R3
X1=X0*C(1,1)+Y0*C(1,2):Y1=X0*C(2,1)+Y0*C(2,2):Z1=X0*C(3,1)+Y0*C(3,2)
G7=T*G1+G2:G7=(G7-INT(G7))*P2:S7=-SIN(G7):C7=COS(G7)
X=+(X1*C7)-(Y1*S7):Y=+(X1*S7)+(Y1*C7):Z=Z1:RETURN

parameters:
X5=(X-X9):Y5=(Y-Y9):Z5=Z-Z9:R5=SQR(X5^2+Y5^2+Z5^2)
IF T6<>T THEN R8=((R6-R5)/(T6-T))/86400 
IF T6=T THEN R8 = -9e+09
R6=R5:T6=T
Z8=+(X5*C8*C9)+(Y5*S8*C9)+(Z5*C9)
X8=-(X5*C8*S9)-(Y5*S8*S9)+(Z5*C9):Y8=(Y5*C8)-(X5*S8)
S5=Z8/R5:C5=SQR(1-S5*S5):E9=ATN(S5/C5)/P0
A9=ATN(Y8/X8)/P0
W5=(ATN(Y/X))*57.3
IF X<0 THEN W5=180-W5
IF X>0 AND Y<0 THEN W5=-W5
IF  X>0 AND Y>0 THEN W5=360-W5
IF X=0 AND Y>=0 THEN W5=270
IF X=0 AND Y<0 THEN W5=90
IF X8<0 THEN A9=A9+180
IF X8>0 AND Y8<0 THEN A9=360+A9
IF X8=0 AND Y8>=0 THEN A9=90
IF X8=0 AND Y8<0 THEN A9=270
B5=Z/R:L5=ATN(B5/(SQR(1-B5^2)))/P0
RETURN